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The authors in their previous papers obtained compact, arbitrarily accurate expressions for two- 
center one- and two-electron relativistic molecular integrals expressed over Slater-type orbitals. In 
this present study, the accuracy limits of given expressions is examined for three-center nuclear 
attraction integrals, which are the first integral set do not have analytically closed form relations. 

They are expressed through new molecular auxiliary functions obtained via Neumann expansion 
of Coulomb interaction. The numerical global adaptive method is used to evaluate these integrals 
for arbitrarily values of orbital parameters, quantum numbers. Several methods, such as Laplace 
expansion of Coulomb interaction, single-center expansion, Fourier transformation method, have 
been performed in order to evaluate these integrals considering the values of principal quantum 
numbers in the set of positive integer numbers. This is the first attempts to study the three-center 
integrals without any restrictions on quantum numbers and in all ranges of orbital parameters. 
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I. INTRODUCTION 

The LCAO-SCF [1] method is generally employed for 
molecules, in which molecular wave functions taken to 
be linear combinations of atomic basis functions whose 
should possess the cusps condition at the nuclei [2] and 
decay exponentially for large distances [3]. This approach 
leads to use, namely, Slater-type orbitals [4, 5], 

Xnin.{C.r)= ( 1 ) 

\/r(2n-bl) 

here, Yim are complex or real spherical harmonics iY*^ = 
Yi-rn]Yim = Sim) differs from the Condon—Shortley 
phases by sign factor (—1)™ [6-8], r(z) are gamma func¬ 
tions [9], {n,l,m} are the principal, orbital, magnetic 
quantum numbers with, n £ K.’*', 0 < I < [nj, —I < 
m < I and [nJ stands for the integer part of n, respec¬ 
tively, in one— and two—electron multi—center molec¬ 
ular integrals. These integrals needs to be calculated in 
spectroscopic accuracy in order to meaningful discussions 
on basis-set expansion methods, Born-Oppenheimer en¬ 
ergy, vibrational frequency calculations. The difficulty of 
finding analytically closed form relations, however, for 
molecular integrals have more than two-center referred to 
as The bottleneck of quantum chemistry [10], have been 
greatest obstacle since Slater-type orbitals have no simple 
addition theorem; relations for products of two Slater- 
type orbitals centered on different positions not available 
in compact form [11]. 

The Slater-type orbitals are obtained by simplification 
of Laguerre functions in hydrogen—like orbitals [12] by 
keeping only the term of the highest power of r, for in¬ 
teger values of principal quantum number n (ISTOs), 
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where n £ N’*', r(2n -|- 1) = (2n)! and it has been proved 
that they provide extra flexibility for closer variational 
description of atoms and molecules by considering the 
values of n in more general set of number, namely posi¬ 
tive real numbers (NSTOs), where n £ R"*". The studies 
on the evaluation of molecular integrals, thus, are per¬ 
formed in two main group: those restrict the principal 
quantum number with integer values, which are practi¬ 
cally used in nonrelativistic molecular electronic struc¬ 
ture calculations [13, 14] and those free them from any 
specification but also reduce the area of applications only 
to investigation of atoms [15-20]. 

The multi—center molecular integrals over ISTOs can 
be evaluated by expansion of Slater-type orbitals through 
complete orthonormal basis functions to a new origin [21- 
25] (see also references therein), 

Xnlm (Cl 

= J™ X ^nTm.n'Z'm'(C,-RAB)Xnq'm'(CUB). (2) 

Ne^OO ^^ ’ 

n'l'm' 

or by expressing them as a finite linear combination of 
B functions through Fourier transform [26-31]. However, 
infinite series representation formulas arising in expan¬ 
sion method require increasing upper limit of summation 
as much as possible to converge to exact values with suf¬ 
ficient decimals (the choice adopted as threshold for the 
total energy in nonrelativistic variational energy calcu¬ 
lation is of order E—03 atomic units, therefore, consti¬ 
tute matrix elements should be accurate to E—10 atomic 
units) and presence of spherical Bessel functions brings 
computational difficulties in Fourier transform method 
since they provoke an oscillation [32-35]. 

The problem of multi—center integrals evaluation by 
the use of NSTOs even much more through insurmount¬ 
able. The Slater type orbitals with noninteger principal 
quantum numbers do not have infinite series represen¬ 
tation formulas; they can not be expanded via complete 
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FIG. 1: Depiction of the coordinates for motion of an 
electron in the field of three stationary Coulomb 
centers, namely A, B, C, where A = {Za,Ya,Xa}, 
B ^ {Zb,Yb,Xb}, C = {Zc,Yc,Xc}, {Z,Y,X} are 
the axes of Cartesian coordinates. 


a6—initio and density functional theory. They are central 
to the understanding of multi-center integrals. They have 
been commonly studied with methods presented above. 
They can also be evaluated through Neumann [40-43] 
and Laplace expansion [44-46] of Coulomb interaction in 
prolate spheroidal coordinates. 

In a new approach the two—center integrals have been 
calculated by the authors for arbitrary values of pa¬ 
rameters and quantum numbers via numerical integra¬ 
tion techniques [47, 48]. The new relativistic molecu¬ 
lar auxiliary functions in prolate spheroidal coordinates 
are presented [48]. They are used to obtain compact 
form relations for two-electron integrals. Afterwards this 
idea adapted to calculate overlap integrals via Fourier 
transform formulas [49]. The same accuracy, 36-digits, is 
achieved in both methods. These are so far only known 
precise calculations for molecular integrals over NSTOs. 
Hence, they are used in this paper to produce benchmark 
values for three-center one-electron molecular Dirac inte¬ 
grals as a first time in the literature. The Mathematica 
programming language [50] is utilized for both analytical 
and numerical calculations. 


II. THREE-CENTER NUCLEAR ATTRACTION 
INTEGRALS 


orthonormal basis functions since power series for a func¬ 
tion such as zP, z £ C and p G R/Nq are not ana¬ 
lytic at the origin [36, 37], where the symbols C, R, 
No used to denote the sets of complex, real and natu¬ 
ral numbers, respectively. It should be noted that, this 
also eliminates possibility of applying binomial expan¬ 
sion theorem in order to evaluate the two-center integrals, 
those are analytically closed form relations may obtain. 
Therefore, in mathematical point of view evaluation of 
multi—center molecular integrals using noninteger prin¬ 
cipal quantum numbers in Slater-type orbitals is an open 
question. It is far more than better representation of elec¬ 
tronic wave—function in nonrelativistic electronic struc¬ 
ture calculations it is also directly related with solution 
of the Dirac equation in algebraic approximation. The 
basis functions to be used in solution of matrix form of 
the Dirac equation are obtained analogously to L-spinors 
[38, 39] which are related to the Dirac hydrogenic solu¬ 
tions. Their explicit form include power functions r''', 

with, Z is nuclear charge, c is speed of light, k = 
±1, ±2, ±2,..., respectively. They can only be represent 
by finite summation of Slater-type orbitals with noninte¬ 
ger principal quantum numbers. 

In particular, the three—center integrals are the first 
set of multi—center integrals do not have analytically 
closed form relations. They have a fundamental im¬ 
portance in the study of molecular systems through 


Taking into account Fig. 1, where depiction of coordi¬ 
nates are given for one electron in a triangular confor¬ 
mation, the three-center nuclear attraction integrals are 
defined as follows, 

^nlm^n'l'm' (CjC^Rab^Rac) 

X*uim \ - rXn'Z'm' U' ,r - Rab] dV, 


with. A, B, C are three arbitrary points of the euclidian 
space, Rab = AB, Rac = AC. 

The Neumann expansion for l/jr-in prolate 
spheroidal coordinates (^, (/>), where l<5<oo, —1< 

:/ < 1, 0 < (/) < 27r [51], 

1 1-.M (^ ~ l-^D- 

\r-RAc\ lm 

xPr(e<)Qf'(?>) 

xYt‘{vcAc)Yt{c<l^r, (5) 

here, are first and second kind associ¬ 

ated Legendre functions [9], {C<)'f>} refers to lesser and 
greater of respectively, is utilized in order to ob¬ 

tain expressions for the three-center nuclear attraction 
integrals over Slater-type orbitals, where the principal 
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quantum numbers are free from specifications; 

l'm' {C,C,Rab,Rac) 

_ 4v^ M 

X yf (i/C,</>c) {sf 

+ , (6) 

here, 

Mnn'{<:X\R) = 

+1/2 

[r( 2 n + l)r( 2 n' + V 2 / 

are the normalization constants and coefiicients [52] 

1 / , \ 1/2 
^^;.+ =7^(2-|i7“|) <5m, 


r-LM 


j-OC i^ + l 

{cx',rab,^c)= / 

Jfc J -1 


nlm.n'l'm 


'ic ■1-1 

X e-«[K?+';')'f*-^i5]->^[i(C-C')-RAB] 


— /l + +\— /1- + 

X P/™ —^ iP/-„ ' ^ 


\e\m—m'\ 


1 


I JC /'0\ 

+ 'lM,e|m+m'|i (,oj 


xQf (e)^LM+)de+ (13) 


where, Pi|m+) are the normalized associated Legendre 
functions. 

The auxiliary functions ,K.^^ can be defined in 
a simpler form through relation given for product of 
two normalized associated Legendre functions centered 
on points A, B in prolate spheroidal coordinates [40], 

A_ X + A^ 

- 1)(1 - +)] Vixicos eA)Vl'X'{cos 9 b) 

I l' ct+/3—2A.—A—A 

= y If y sifiixi'>.’A) 

a=-(2A-A) /3=X' q=0 

(+)^ 


(^ + - j,)/3 


(14) 


are the integration over azimuthal angle. The symbol e 
may have the value ±1 and is determined by the product 
of the signs m and m' (the sign of zero is regarded as 

I / 

positive). The symbols may have the values ±1 

and 0: if among the indices m, m' and m ± m' there 
occurs a value equal to zero, then is also zero; if 

all the indices differ from zero, rjm^ = ±1 and the sign 
is determined by product of the signs m, m', m ± m'. 
Thus, the coefficients A^ differ from zero only with the 
values \M\ = jm — m'j, \M\ = jm + m'j. 

Since it is assumed axes of prolate spheroidal coordinate 
system centered on A, B substitutions in Eqs. (4, 5) can 
be written as follows [53], 


which is obtained in explicit form by taking advantage of 
binomial expansion theorem, 

(x + a)^^{x — 

N 1 +N 2 

= Fs{NuN2)x^^+^^-^a\ (15) 


s=0 


with, 


cos0yi = 


1 + 


COSC/r = 


l-^B 



^ rA+rB , 
Rab ’ 

rA - rs 
^ , 

Rab 

(9) 

by. 

XA = 

^AB /> I A 

2 + 

Rab / * \ 

rB= 2 A v), 

(10) 


+ = 

Kac+Rbc 

Rac—Rbc 

(11) 


Rab 

Bn = -■ 

Rab 



The coefficients occurring in Eq. (14) are determined 
9lp{l\ ^'A) = glpilX, l'\)F,{a + A, /3 - A) (16) 

A 

g%{l\l'\) = Y.^-iyPs{\)D^^+ 2 X- 2 sDl\ (17) 

s=0 


The integrals are the auxiliary functions and 

they are defined as, 




2l + lFi{l + \) 


1/2 


rrLM 
nlm.n' I'm‘ 


(C, c', RAB^ic) = f\i + l^)^{i-l^)^' 


2 Fx{l) \ 

^F^^i_py2{l)Fp-x{l + P), (18) 


X Vim ( ^ ^ 1 Vl' 




1- + 
i-u 


xP^^^{OVLM{y)d^du ( 12 ) 


where, A = [mj, A' = Jm'] and the quantities Fs{N,N') 
are the generalized binomial coefficients. They are given 
as, 

FyN,N') = Y,{-y"'Fs-s'{N)F,,{N') (19) 
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with, i [(s — -/V) + |s — A^|] < s' < min{s, N) and Fs{N) 
are binomial coefficients indexed by N and s is usually 

written ^ ^ ^ j respectively. 

The Eqs. (12, 13) are, therefore, obtained as follows. 


(C, C', RaB,^c) 

-rLM.: 


= E5a/3aA,rA';A) 

aPq 


iC, C,Rab,^c) 
{C,C\RaB.^c) 


with, 


( 20 ) 


Jn-a,u'-l3 (C, C', RaB , ) 

.o] 




p|A| 

Qf (0] 


'PL\A\{y)d^dv. (21) 


There are no known convergent series representation for¬ 
mulas, free from specifications on parameters for power 
functions such as (^ + i^)^b (C — *^)^^, {-^ 1 ,-^ 2 } € K, 
yet and that poses an obstacle to analytically reduce 
the auxiliary functions to one variable 

auxiliary functions introduced in [42]. Thus, 
the solution should be obtained on the basis of numer¬ 
ical methods. Note that, taking advantage of binomial 
expansion method for terms containing the angular part 
of Slater-type orbitals in order to simplify the expres¬ 
sions increases the number of integrals should be numeri¬ 
cally calculated. In Mathematica programming language 
instead of using Eq. (20), direct computation of Eqs. (12, 
13) are faster. The given relations for auxiliary functions 
in Eq. (21) are calculated by using different expressions 
of Legendre polynomials and compared according to com¬ 
putational time in Fig. (2) as sample. The discussions on 
results are made in the next section. 


rithm described in [47] has been incorporated into a com¬ 
puter program written in the Mathematica programming 
language with the included numerical computation pack¬ 
ages for solving Eqs. (6, 12, 13). The Mathematica pro¬ 
gramming language can handle approximate real num¬ 
bers with any number of digits and it is suitable for 
benchmark evaluation. It is also provides a uniquely in¬ 
tegrated and automated environment for parallel com¬ 
puting. It is allows us to compute the formulas including 
summations using all cores of PC effectively via Parallel- 
Sum command instead of Sum. Note that, in this study 
all results are given in atomic units (a.u.). The calcula- 


TABLE I: The values for auxiliary functions defined in 
Eq.(21) with pi = ^ {( + (') R^b, P 2 = ^ {C - C) Rab 
and Ni,N 2 € K"*". 


L A q Ni Ni Pi p2 
0 0 0 1.0 1.0 1.5 1.5 


1 0 2 3.0 2.0 4.5 0.5 


2 2 2 3.0 2.0 4.5 0.5 


3 2 5 9.0 4.0 22.5 0.1 


5 4 6 15.0 3.0 4.0 0.1 


0 0 0 1.2 1.5 1.5 1.5 


1 0 2 3.3 2.4 4.5 0.5 


2 2 2 3.5 2.5 4.5 0.5 


3 2 5 9.9 4.1 22.5 0.1 


5 4 6 15.2 3.6 4.0 0.1 


5 4 6 15.2 3.6 0.1 4.0 


5 4 6 15.2 3.6 0.1 4.0 


Results 

3.62319 79582 17897 45490 E-01 
1.75859 65139 47296 72718 E-01 
1.02307 57195 13525 65648 E-03 
1.20185 35031 71549 79713 E-03 
1.26482 46553 60927 73809 E-02 
3.05572 36905 83528 19812 E-04 
1.54156 95966 91532 03328 E-11 
9.77333 04203 13413 96225 E-18 
7.45864 92397 67729 51682 E-h06 
1.65225 27526 05586 45258 E-h05 
4.93516 66112 08595 80377 E-01 
3.72444 04752 38238 19870 E-01 
7.02522 66592 59862 42272 E-04 
4.07093 04267 92761 66995 E-05 
2.01392 20090 23026 29191 E-02 
6.86518 34228 55977 92648 E-04 
2.74437 86677 61624 01320 E-11 
2.53854 67628 59329 46820 E-17 
9.75947 26622 39909 40431 E-h06 
4.76137 54661 40175 09867 E-h05 
9.01557 46105 79752 64049 E-h08 
9.49361 57667 79094 34667 E-h37 
6.67509 82662 40505 33704 E-h08 
1.34769 33998 41265 16644 E-h36 


III. RESULTS AND DISCUSSIONS 

The literature currently, lack of benchmark values 
for multi—center integrals when Slater-type orbitals 
are used. Recently, a robust numerical Global—adaptive 
strategy with Gauss—Kronrod extension has been ap¬ 
plied for two—center integrals through prolate spheroidal 
coordinates and fourier transform method in [47-49]. 
Benchmark results have been presented for them. The 
hermitian properties are thus, represented correctly free 
from specification on quantum numbers, orbital param¬ 
eters and internuclear distances. In this study it is ex¬ 
tended for solution of three—center integrals. The algo- 


tion results are presented in Tables I—VII and Fig.2 for 
arbitrary values of quantum numbers, orbital parameters 
and internuclear distances. The comparisons are made 
with expansion methods which are given for expansion 
of wave function in Eq. (2) and for charge density ex¬ 
pansion to same center by following formula [54-56] 


Pnlm,n'Cm' iCi ^^ C 
l+l' I" 

— ^ ^ ^ ^ n'I'm\n+n*— 11"m" iCj C t 

= I m" =—l" 

^ Xn+n'— ( 22 ) 
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They are useful to reduce the three—center integrals to 
basic nuclear attraction integrals, 

-.If 1 

J^Xr{z,RBc) = / X*K\r{^,rB) — dvi- (23) 

V47r J rc 

The Eqs. (2, 22) are used to transform the wave function 
centered at A to a wave function centered at B then, to 
transform the charge density centered on same positions 
to a single wave function, respectively. Here, z = C, + 
Cl Rbc = The resulting basic nuclear attraction 
integrals are calculated by the following formula, [57], 

JnXri^Z^R^ = 

2^ f2r{K + X + 2) 1 

2X + l\l z x/T{2k + 1) {zR)^+^ 

T{k + X + 2,zR) {zRC^+^r{K-X + l,zR)\ 

r(K +A + 2) r(K +A + 2) J 

xYxAO,^), (24) 

where, r(n,m) are incomplete gamma functions [9]. 

In Fig. 2 the Eq. (21) is investigated according to com¬ 
putational time in Mathematica programming language. 
Mathematica includes all the common special functions 
of mathematical physics. It also provides easy way of 
computing them precisely. Here, explicit formula (EE) 
[58], 

Vixix) = {l-x^yYl (25) 

k 



— — 

OlX — -■ 


21 + 1 


21 [2Fx{l)Fx{l + X)\ 

X (-l)'=Ffe(A + k)Fi_C‘2l - 2k)Fi_x-2k{l - k), (26) 


where, 0 < k < E , recurrence relation formula 
(RE) [59] of Legendre polynomials are compared with 
Mathematica built-function (ME) LegendreP[n,m,x]. 
They are presented with red, blue, green lines in Fig. 2, 
respectively. It can be seen from this figure, the direct 
use of Mathematica buit-function, given for computing 
of Legendre polynomials, in numerical integration of Eq. 
(21) provide the results faster then explicit or recurrence 
relation formulas. 

The results for calculation of Eq. (21) are also presented 
in Table I with integer and noninteger values of prin¬ 
cipal quantum numbers. The first, second rows are ob¬ 
tained from calculation and functions, re¬ 
spectively. The auxiliary functions defined in 

[42] for three-center integrals are special case of Eq. (21). 
Hence, we believe an importance of present the results for 
general form of Eq. (21). 

The results for Eqs. (12, 13) and Eq. (20) are presented 
in Tables H, HI. They are given in first, second rows 
for and auxiliary functions, re¬ 

spectively. Note that, the Eq. (20) is only differ from 


Time(s) 



FIG. 2: CPU time for computation of auxiliary 

function in Eq.(21) according to methods used for 
calculation of the Legendre polynomials, where, 
Mathematica built-function (ME), recurrence relation 
formula (RE), explicit formula (EE) and, L = 3, A = 1, 
g = 0,Ni = 3,N2 = 2, Pi = 2.5, p2 = 1.5. 


Eqs. (12, 13) in that, the normalized associated Legen¬ 
dre polynomials on right-hand side are expanded via Eq. 
(14) and the numerical global adaptive method is per¬ 
formed to remaining parts. Performing the calculations 
in Mathematica programming language for such formulas 
contain summations is disadvantageous in terms of cal¬ 
culation time. The results in Tables H, HI shows that, 
numerical Global adaptive method with Gauss-Kronrod 
extension can be used for computation of Eqs. (12, 13) 
which is eliminate necessity applying binomial expansion 
theorem. 

In Tables IV, V, VI the results obtained for three- 
center integrals are presented for upper limit of summa¬ 
tion L is L = 30. The correct digits are underlined. The 
digits in bold indicate the convergence property of used 
method. In first and second rows benchmark results ob¬ 
tained from numerical global adaptive method and re¬ 
sults those found in the literature are presented, respec¬ 
tively. Later rows are the results obtained from expansion 
of the wave-function method and they are in complete 
agreement with ones from E. §ahin (personal commu¬ 
nication), where the upper imit of summations Nf, are 
given in parenthesis. The expansion of the wave-function 
method is tested up to upper limit of summation Vg, 
Ne = 160. It is observed that, the results hardly con¬ 
vergent with 10—digits for given quantum numbers and 
orbital parameters in table V. In other tables the con¬ 
vergence remains between 5—digits and 10—digits. Note 
that, it is necessary to take into account eight summa¬ 
tion and four of them should be infinite in expansion of 
the wave-function method if NSTOs are used. The values 
presented in Tables IV, V, VI for STOs clearly demon¬ 
strate pointlessness of such an attempt. 

On the other hand in our previous papers [47, 48] it have 
been proved that the numerical Global adaptive method 
with Gauss-Kronrod extension is able to give benchmark 
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values. In particular for Table VII the results are pre¬ 
sented for different upper limit of summation L appears 
in Eq. (6). They differs from upper limit of summation 
Ne used in expansion of STOs in that they are given in 
brackets. Convergence property of Eq. (6) is examined 
in this table. It is found that, by increasing the upper 
limit of summation the results are convergence to exact 
values. The results obtained up to upper limit of summa¬ 
tion L is L = 40. It is achieved to 25—digits accuracy by 
determining the upper limit of summation L is L = 30 
accordingly, there is no necessity of performing calcula¬ 
tions with upper limit of summation higher than L = 30 
unless more precise results required for a given values 
of parameters. It should be point out that, the summa¬ 
tion appears in Eq. (6) should not be regarded as having 
same characteristic with summation arising in expansion 


of STOs. It is based on expansion of spherical harmon¬ 
ics which have form a complete set of orthonormal func¬ 
tions. Any square-integrable function can be expanded 
as a linear combination of spherical harmonics. The con¬ 
vergence problems arising in expansion NSTOs can not 
exist in our method. Without any computational diffi¬ 
culty by increasing the upper limit of summation L can 
be achieved to desired accuracy rapidly. 
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TABLE II 

; The values for auxiliary functions defined in Eqs. (12,13) and Eq. (20) with = 1 (^ + (') Rab, 

P2 = 5 (C — CO ^ab and n, n' G N+. 

L 

M 

n 

1 

m 

n' 

V 

m! 

Pi 

P2 

Eqs. (12,13) and Eq. (20) 

0 

0 

1.0 

0 

0 

1.0 

0 

0 

1.5 

2.5 

2.90399 07988 62162 00938 68706 03684 E-01 
1.47689 95097 54804 95414 54214 62847 E-01 

2 

0 

3.0 

0 

0 

2.0 

0 

0 

0.01 

2.5 

3.42865 83511 73345 68718 16357 54166 E-00 
8.97033 20343 34376 52929 54799 07851 E+05 

5 

0 

9.0 

0 

0 

4.0 

0 

0 

22.5 

0.1 

6.61192 17510 24844 44964 35310 16794 E-12 
1.27733 49592 55934 35439 28200 66327 E-23 

10 

0 

12.0 

0 

0 

8.0 

0 

0 

40.0 

0.001 

1.73471 85107 61691 61997 46481 63897 E-20 
1.87470 00384 84104 39018 32142 95122 E-41 

1 

0 

3.0 

1 

0 

2.0 

1 

0 

2.3 

4.5 

1.06208 64400 83910 27727 15664 09620 E-00 
2.45315 16908 16648 21895 84639 84955 E-01 

6 

0 

6.0 

4 

0 

5.0 

3 

0 

0.001 

0.001 

2.25821 42942 16571 75095 60318 57985 E+46 
2.48928 81670 28192 48883 41015 71291 E+56 

1 

0 

2.0 

1 

1 

2.0 

1 

1 

0.01 

0.01 

1.39224 44783 13768 61934 60782 72901 E-02 
1.63293 63668 41775 85565 65259 43851 E+03 

3 

0 

6.0 

4 

3 

5.0 

3 

3 

8.5 

0.1 

1.22009 58147 71182 05449 05454 44137 E-04 
1.63709 36934 93917 83173 07173 88885 E-08 

3 

2 

5.0 

3 

3 

3.0 

2 

1 

1.5 

1.5 

8.36879 35046 74242 89056 27574 95673 E+01 
2.18611 99801 01738 31430 63739 74869 E-00 

3 

2 

10.0 

4 

4 

8.0 

2 

2 

0.1 

9.0 

8.73612 92750 76127 55419 06451 97990 E+06 
1.89278 99267 51317 93832 38020 79768 E+06 







TABLE III; The values for auxiliary functions defined in Eqs. (12,13) and Eq. (20) with = 1 (^ + (') Rab, 

P2 = 5 (C — CO Rab and n, n' G IR+. 
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L M n I m n' I' m' Pi p 2 

0 0 1.1 0 0 1.5 0 0 1.5 1.5 


1 0 1.3 0 0 1.8 0 0 0.5 0.1 


2 0 3.2 0 0 2.4 0 0 0.01 2.5 


5 0 9.5 0 0 4.8 0 0 22.5 0.1 


10 0 12.5 0 0 8.9 0 0 40.0 0.001 


1 0 3.3 1 0 2.6 1 0 2.3 4.5 


6 0 5.9 4 0 6.1 3 0 0.001 0.001 


1 0 2.3 1 1 2.1 1 1 0.01 0.01 


3 0 5.9 4 3 5.1 3 3 8.5 0.1 


3 0 8.8 4 3 10.3 3 3 0.01 0.01 


3 2 4.5 3 3 3.5 2 1 1.5 1.5 


3 2 10.5 4 4 8.1 2 2 0.1 9.0 


Eq.(16,17) and Eq.(19) 

4.09565 63231 31448 87039 55147 40491 E-01 
2.93614 23789 04133 04860 14608 03840 E-01 
3.33243 32418 64279 67130 02542 67866 E-01 
1.02388 22819 29861 91096 69560 06751 E-02 
5.49415 76314 41209 10218 23695 48480 E-00 
7.89858 88444 51044 77570 43625 77546 E+06 
8.85039 89771 63772 65903 44204 90792 E-12 
9.74443 14137 20230 92395 49657 75617 E-23 
1.17276 22934 30019 21937 27274 20564 E-20 
5.51171 88865 58272 19013 29428 46796 E-41 
1.87161 11340 32480 14502 70364 65188 E-00 
6.65100 22118 49834 63694 23667 20358 E-01 
7.38511 70805 47763 29724 30676 68341 E+46 
1.24500 66118 34410 00379 04215 62074 E+60 
1.85417 21104 44330 08693 85023 65883 E-01 
1.40783 85008 42084 72060 75228 55754 E+04 
1.13313 80933 64382 91558 05780 65710 E-04 
1.56774 69588 19924 39645 11990 75063 E-08 
4.97973 96013 63878 34830 92006 86281 E+04 
7.95070 98832 40529 88293 62402 30270 E+42 
1.12197 95987 12589 22596 64533 39225 E+02 
2.51548 44999 37750 08034 41787 62496 E-00 
1.10758 39108 61744 40668 98531 33050 E+07 
1.32690 39575 82905 48879 15023 14344 E+29 


TABLE IV: 

The values for three-center 

integrals, where position of the nuleus A, B, C in cartesian coordinates 




{X,Y,Z}: 

A={0,0,0}, B={0,0, 

-2.0143}, C={0,0, 

—4.1934}, respectively. 

n 

1 

m 

C 

7l' 

1' 

vn! 

C' 

Results 









2.94549 60536 73751 14101 41604 E-02 









2.94549 6054 E-02 

2.92367 19340 54421 10253 45246 E-02 (5) 

1.0 

0 

0 

1.24 

1.0 

0 

0 

5.67 

2.95420 14999 28792 85854 63190 E-02 (10) 
2.94569 30497 61154 50027 23145 E-02 (20) 
2.94549 53848 92871 61006 54907 E-02 f40') 
2.94549 58796 89501 79906 76961 E-02 f80i 
2.94549 60532 13114 00159 34423 E-02 (160) 

1.1 

0 

0 

1.24 

1.1 

0 

0 

5.67 

3.34089 64668 60903 49505 70936 E-02 

1.2 

0 

0 

1.24 

1.2 

0 

0 

5.67 

3.74811 24252 41741 71875 67668 E-02 

1.3 

0 

0 

1.24 

1.3 

0 

0 

5.67 

4.16291 70512 36793 53772 11297 E-02 









1.60664 60408 09377 56494 10579 E-01 

1.60664 6078 E-Or 

1.57770 82722 80008 33479 04172 E-01 (5) 

1.0 

0 

0 

1.24 

2.0 

0 

0 

1.61 

1.60084 83584 29304 06752 43872 E-01 (10) 
1.60753 33306 97032 54324 80156 E-01 (20) 
1.60673 68838 68511 26148 86494 E-01 (40) 
1.60659 82475 67392 07711 26316 E-01 (80) 
1.60664 10673 25488 09376 70768 E-01 (160) 

1.1 

0 

0 

1.24 

2.1 

0 

0 

1.61 

1.69763 82691 19573 36045 53004 E-01 

1.2 

0 

0 

1.24 

2.2 

0 

0 

1.61 

1.78213 08867 60019 08900 15372 E-01 

1.3 

0 

0 

1.24 

2.3 

0 

0 

1.61 

1.86000 31133 33942 68781 62339 E-01 


[ 11 ] 
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TABLE V: The values for three-center integrals, where ( = (' = 2.0; position of the nuleus A, B, C in cartesian 
coordinates {X,Y,Z}: A={0,0, 0}, B={0,0, 6}, C={0,0, —7}, respectively. 


n 

1 

m 

n' 

v 

m' 

Results 


4.53377 50011 42666 45050 53528 E-04 
4.53377 50011 426 5 E-04^ 

4.53377 50011 15138 130 E-04'" 


4.49978 10366 11814 81692 67046 E-04 (5) 


2.0 

0 

0 

2.0 

0 

0 

4.53371 93448 02014 78421 78747 E-04 (10) 
4.53377 50006 71569 59704 90361 E-04 (20) 
4.53377 50013 47449 39243 55321 E-04 (40) 
4.53377 50011 98380 66714 63012 E-04 (80) 
4.53377 50011 56735 48522 57686 E-04 (160) 

2.1 

0 

0 

2.1 

0 

0 

5.51824 66271 98490 20314 29459 E-04 

2.2 

0 

0 

2.2 

0 

0 

6.66626 60479 98863 93114 16105 E-04 

2.3 

0 

0 

2.3 

0 

0 

7.99547 53772 28315 43332 78068 E-04 

2.0 

1 

0 

2.0 

0 

0 

6.49989 25189 10135 10589 99466 E-04 
6.49989 25189 1011 E-04"- 

2.1 

1 

0 

2.1 

0 

0 

7.90609 20760 42277 16143 10272 E-04 

2.2 

1 

0 

2.2 

0 

0 

9.54293 49790 09494 09113 15449 E-04 

2.3 

1 

0 

2.3 

0 

0 

1.14342 80666 86023 44617 17750 E-03 

2.0 

1 

0 

2.0 

1 

0 

9.69666 68121 64802 10682 05157 E-04 
9.69666 68121 6477 E-04" 

2.1 

1 

0 

2.1 

1 

0 

1.17838 66960 60512 32340 53458 E-04 

2.2 

1 

0 

2.2 

1 

0 

1.42089 42554 34926 39305 88294 E-04 

2.3 

1 

0 

2.3 

1 

0 

1.70054 57496 59994 51520 55390 E-03 


[41] 

[43] 


TABLE VI: The values for three-center integrals, where C = C' = 1.6; position of the nuleus A, B, C in cartesian 
coordinates {X,Y,Z}: A={0,0, 0}, B={0,0, 3}, C={0,0, —3}, respectively. 


n 

1 

m 

n' 

i' 

m' 

Results 


7.41579 46662 21323 37855 75053 E-02 
7.41579 46662 2133 E-02^ 

7.41579 46662 21323 37 E-02'" 


7.39466 24791 10067 86601 96990 E-02 (5) 


2.0 

0 

0 

2.0 

0 

0 

7.41595 39598 23811 62963 42293 E-02 (10) 
7.41574 64572 83043 94945 93805 E-02 (20) 

7.41578 73240 73721 74843 59312 E-02 (40) 

7.41579 22349 71134 02044 81338 E-02 (80) 
7.41579 51272 77100 02097 74777 E-02 (160) 

2.1 

0 

0 

2.1 

0 

0 

7.92564 04306 46454 51287 28675 E-02 

2.2 

0 

0 

2.2 

0 

0 

8.42445 84406 53803 96918 57249 E-02 

2.3 

0 

0 

2.3 

0 

0 

8.90946 29369 56705 18394 86907 E-02 

2.0 

1 

0 

2.0 

0 

0 

6.77544 93679 78440 06490 13878 E-02 
6.77544 93679 7845 E-02" 

2.1 

1 

0 

2.1 

0 

0 

7.10247 82121 12668 99509 25056 E-02 

2.2 

1 

0 

2.2 

0 

0 

7.39664 78499 31562 91249 65250 E-02 

2.3 

1 

0 

2.3 

0 

0 

7.65548 51830 78592 48565 84851 E-02 

2.0 

1 

0 

2.0 

1 

0 

6.75368 18945 49233 36334 86913 E-02 
6.75368 18945 4926 E-02" 

2.1 

1 

0 

2.1 

1 

0 

6.91775 52123 72612 96334 96130 E-02 

2.2 

1 

0 

2.2 

1 

0 

7.02466 73478 35283 82859 84242 E-02 

2.3 

1 

0 

2.3 

1 

0 

7.07260 26896 46234 08137 55301 E-02 


[41] 
*" [43] 





















11 


TABLE VII: The values for three-center integrals, where C = 3.6, C' = 1.6; position of the nuleus A, B, C in 
cartesian coordinates {X,V,Z}: A={0,0, 0}, B={0,0, 3}, C={3,0, 3}, respectively. 


n 

1 

m 

r\! 

1' 

m! 

Results 

2.0 

0 

0 

2.0 

0 

0 

2.72309 91701 60088 20662 71607 E-02 [1] 
2.70276 01923 98151 89286 78890 E-02 [51 
2.70272 90219 79709 38269 45472 E-02 [101 
2.70272 90187 20151 76256 32122 E-02 [111 
2.70272 90189 88329 48522 85282 E-02 [121 
2.70272 90190 09862 41008 43703 E-02 [151 


2.70272 90190 09906 8540 6 45442 E-02 [20] 
2.70272 90190 09906 85401 14157 E-02 [30] 
2.70272 90190 09906 85401 14157 E-02 [40] 
2.70272 901 89 8 E-02‘" 

2.70272 90190 2 E-02'" 


2.1 

0 

0 

2.1 

0 

0 

2.99614 53594 03007 59134 38353 E-02 [101 
2.99614 53558 58718 15620 66646 E-02 [151 
2.99614 53558 58738 85669 00599 E-02 [301 
2.99614 53558 58738 85669 00599 E-02 [40] 

2.2 

0 

0 

2.2 

0 

0 

3.30034 60724 25497 07263 17130 E-02 [101 
3.30034 60682 94929 39385 11310 E-02 [151 
3.30034 60682 94915 87198 38700 E-02 [301 
3.30034 60682 94915 87198 38698 E-02 [40] 

2.3 

0 

0 

2.3 

0 

0 

3.61354 34150 29555 79331 66393 E-02 [101 
3.61354 34103 21082 87274 34304 E-02 [151 
3.61354 34103 21023 85552 77843 E-02 [301 
3.61354 34103 21023 85552 77839 E-02 [40] 

2.0 

1 

0 

2.0 

0 

0 

2.15396 79080 46322 59289 24910 E-02 [101 
2.15396 78874 84626 90451 70130 E-02 [151 
2.15396 78874 84896 61732 77914 E-02 [301 
2.15396 78874 84896 61732 77914 E-02 [401 
2.15396 78873 E-02"' 

2.1 

1 

0 

2.1 

0 

0 

2.37228 09111 80051 67464 97075 E-02 [101 
2.37228 08887 59565 69293 67888 E-02 [151 
2.37228 08887 59514 88542 81485 E-02 [301 
2.37228 08887 59514 88542 81456 E-02 [40] 

2.2 

1 

0 

2.2 

0 

0 

2.59362 19462 42909 94539 77789 E-02 [101 
2.59362 19225 73538 25384 87321 E-02 [151 
2.59362 19225 73080 19834 68030 E-02 [301 
2.59362 19225 73080 19834 67975 E-02 [40] 

2.3 

1 

0 

2.3 

0 

0 

2.81580 26527 44329 75185 31334 E-02 [101 
2.81580 26286 47036 51507 08634 E-02 [151 
2.81580 26286 46091 20634 90123 E-02 [301 
2.81580 26286 46091 20634 90048 E-02 [40] 

2.0 

1 

0 

2.0 

1 

0 

3.32260 90918 27546 98165 04057 E-02 [101 
3.32260 90468 18302 55100 75237 E-02 [151 
3.32260 90468 19031 39539 62493 E-02 [301 
3.32260 90468 19031 39539 62493 E-02 [401 
3.32260 90465 E-02" 

2.1 

1 

0 

2.1 

1 

0 

3.64394 52557 31327 16654 23236 E-02 [101 
3.64394 52049 85248 40191 32912 E-02 [151 
3.64394 52049 85324 78984 05227 E-02 [301 
3.64394 52049 85324 78984 05166 E-02 [40] 

2.2 

1 

0 

2.2 

1 

0 

3.96639 94790 03807 01164 82993 E-02 [101 
3.96639 94232 71120 86583 47853 E-02 [151 
3.96639 94232 70337 92234 14511 E-02 [301 
3.96639 94232 70337 92234 14397 E-02 [40] 

2.3 

1 

0 

2.3 

1 

0 

4.28635 93143 98419 21348 89833 E-02 [101 
4.28635 92548 29750 51077 50463 E-02 [151 
4.28635 92548 27904 76970 96420 E-02 [301 
4.28635 92548 27904 76970 96264 E-02 [401 


[41] 

[43] 



















